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Abstract. - We present a new procedure that can identify and measure the critical temperature. 
This method is based on the divergence of the relaxation time approaching the critical point in 
quenches from infinite temperature. We introduce a dimensionless quantity that turns out to 
be time-independent at the critical temperature. The procedure does not need equilibration and 
allows for a relatively fast identification of the critical temperature. The method is first tested 
in the ferromagnetic Ising model and in the two dimensional EA model and then applied to the 
one-dimensional Ising spin glass with power-law interactions. Here we always find a finite critical 
temperature also in presence of a uniform external field, in agreement with the mean-field picture 
for the low temperature phase of spin glasses. 



The identification of a critical temperature T c is a fun- 
damental characterization of statistical systems. This, in- 
deed, allows one to construct the phase-diagram of the sys- 
tem and to obtain insights in the underlying relevant phys- 
ical mechanisms. In many cases the existence/absence of 
a phase-transition discriminates among different pictures 
for a given system. A major example is the long-standing 
question on the nature of the spin-glass phase in finite 
dimensional systems. According to the replica symmetry 
breaking scenario [I], a transition line should exist, the 
so-called de Almeida-Thouless (AT) line, separating the 
paramagnetic from the spin glass phase in the tempera- 
ture vs magnetic field phase diagram [T][2] . More precisely, 
the spin glass phase is not destroyed by applying an ex- 
ternal field. Conversely, other theories, as the "droplet 
picture" [3], predict no AT line. Therefore, the identifica- 
tion of a critical temperature in the presence of an external 
perturbation discriminates between the different theories. 

The critical point is usually identified with the tempera- 
ture T — T c where the equilibrium correlation length £(T) 
diverges. This property reflects on the behavior of the 
order parameter correlation function C(r) whose asymp- 
totic decay changes from exponential to algebraic when 



approaching T c . Such a study is often hindered by the 
fact that also the relaxation time t eq diverges for temper- 
atures close to T c and it is only possible to equilibrate 
systems of small size L. However, it is possible to ex- 
trapolate the critical temperature in the L — > oo ther- 
modynamic limit from the behavior of finite systems by 
means of well established methods such as Finite Size Scal- 
ing (FSS) |3]. The most used method for the identifica- 
tion of T c , generally known as Phenomenological Renor- 
malization, consists in the introduction of an appropriate 
dimensionless quantity that is expected to cross exactly 
at the critical point [5]. In the specific case of magnetic 
systems, one usually measures j6H2] a finite size correla- 
tion length £(L,T) at the temperature T, as for instance 

£(L,T) 2 = (tfdr r 2 C(r,T)) / J Q L drC(r,T) [TO]. The di- 

mensionless quantity X(L, T) = £(L, T)/L is then plotted 
over T for different values of L. According to FSS, one 
expects the following relation 



X(L,T) = f[aT)/L] 



(1) 



implying that X(L, T) becomes L independent when £(T) 
diverges, namely for T = T c . T c is therefore given by the 
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temperature where curves for different L intersect. 

On the other hand one can identify T c exploiting the 
intrinsically non-equilibrium nature of the critical point 
as, for instance, implemented in Short Time Dynam- 
ics (STD) [IT], Non-Equilibrium Relaxation (NER) [T2] 
or other methods which use both equilibrium and non- 
equilibrium measurements |13j . These are substantially 
based on the dynamical scaling hypothesis assuming the 
existence at the time t of a typical dynamic length L(t). 
The divergence of the relaxation time at T c then reflects 
on a power law temporal decay of some observablcs such 
as the magnetization. 

In this letter, exploiting dynamical scaling as in NER 
methods, we propose a procedure to identify T c that uses 
phcnomcnological renormalization techniques with non- 
equilibrium time dependent quantities. More precisely, 
we set L — > oo from the beginning, and replace in Eq. ([1]) 
space distances with times introducing the dimcnsionlcss 
quantity 

1 C dt tC(t, T) 
X(r,T) = - J ° > (2) 

t J dtC(t,T) 

where t is the time since the quench at the temperature 
T from a disordered initial condition, r is a fixed time 
and C(t,T) is the correlation function with the initial 
configuration. For instance, for spin systems, C(t,T) = 
((««(t) - (s i (t)))(s i (0) - (s;(0)))), with mtt) the spin in 
the position i at time t. The average is performed over 
initial conditions, non-equilibrium dynamics at the tem- 
perature T, and quenched disorder if present. We al- 
ways take initial conditions corresponding to equilibrium 
at infinite temperature, and therefore (sj(0)) = and 
C(t,T) = (si(t)si(0)). Let us notice that other quantities 
can be equivalently considered in the definition of X(t, T). 
Here we use the two-time correlation function, that is more 
easily obtained in numerical simulations. However one can 
also use the thermo-remanent magnetization, more suit- 
able in experimental settings. Dynamical scaling, then, 
predicts 



X(T,T)=g[t eq (T)/i 



(3) 



where t eq (T) is the relaxation time for quenches at the 
temperature T. The divergence of t eq (T) approaching T c 
implies that X(t,T) does not depend on r when T = T c . 
Therefore, considering different values of r and plotting 
X(t, T) vs T, one identifies T c as the intersection point of 
the different curves. 

The above procedure does not require to equilibrate the 
system, leading to some advantages with respect to the 
FSS method. Indeed, in experiments, the sample sizes 
always fulfill the thermodynamic limit and one cannot 
explore the L dependence. Conversely, one can measure 
X(t,T) for different time interval r and different T. From 
a numerical point of view, in particular for spin-glasses, 
the equilibration of the system is numerically hard to 
achieve and to check. Indeed, equilibration is very time- 
demanding and indirect checks are always necessary in 
order to verify that a true equilibrium state is attained. 



Moreover, in our approach, by means of a single simulation 
up to the final time tf, one can compute X(t, T) for many 
values of r e [0, tf]. Conversely, in the FSS method one 
has just one X(L,T) for each simulation with a system 
of size L. Another difference relies on the possibility of 
using simple spin-spin correlations as C(t,T) also in dis- 
ordered systems, where, generally, the identification of a 
correlation length necessitates the computation of multi- 
spin correlation functions |14j . The presence of disorder, 
indeed, makes £(L, T) invisible to equal time spin-spin cor- 
relation function. Conversely, C(t, T) is strongly affected 
by t eq (T) . A final remark concerns the possibility to better 
control the influence of scaling corrections in the estimated 
T c . Indeed, a pure power law decay C(t, T) ~ t~ e " must be 
observed at T = T c , leading to X(t,T c ) = (l-0 c )/(2-0 c ). 
The exponent 9 is related to the Fisher-Husc exponent A 
via the relation 6 — X/z, where z is the growth exponent. 

The measured X at the intersection point, therefore, 
gives an estimate for C and corrections to scaling should 
be observed as deviations from this power law decay of 
C(t,T). Let us notice, however, that a precise measure- 
ment of T c can be only obtained for very large r. From 
dynamical scaling, evolution up to a finite time tf corre- 
sponds to equilibration up to a size L(tf), and the accu- 
racy in the determination of T c is then of the same order 
of FSS analysis on system up to size L(tf). In the follow- 
ing we present results involving not very large simulation 
times (~ 12ft, of cpu time for each temperature), that, 
however, are sufficient to identify the critical temperature 
with a reasonable accuracy. For each given system and 
each temperature we consider about 1000 independent re- 
alizations. 

Dynamical evolution is obtained via standard Monte 
Carlo simulations and X(T,t) is obtained after the inte- 
gration of C(t,T) with a time-step of single spin update. 
In all cases we always take different sample sizes L in order 
to check that no finite size effects are present. 

Let us begin by checking our method in cases where 
the critical temperature is well known. In particular, 
we start by considering the Ising model with Hamilto- 
nian % = — J2{ij) Jij s i s j arL d ferromagnetic coupling 
Jij = J in two and three dimensions. In these cases 
T c ~ 2.269J and T c ~ 4.5115J are, respectively, ana- 
lytically and numerically known |15j . In the following we 
always take J = 1 for simplicity. The behavior of C (t, T) 
can be obtained from general arguments [TTJ[T6l[T7] , giv- 
ing C(t,T) - t-9e e -*/t e ,(r) for T > Tcj where q c can 

be related to static and dynamic critical exponents, and 
C(t,T) - t~ e with < 9 C for T < T c . Then one has, for 
large r and 8 C < 1 



X(t,T) 




te q {T)/l 



h)/(2-6c) 
(l-0)/(2-0) 



for T > T c 

for T = T c (4) 

for T < T c . 



More precisely, for T < T c , the dynamics is initially at- 
tracted by the critical point at T c [18] and then converges, 
for large r, to Xq > X(t,T c ). This implies that X(r,i) 



p-2 



Identification of the critical temperature from non-equilibrium time-dependent quantities 





Fig. 1: (Color online). X(t,T) and X2(t,T) are plotted versus 
T for different values of r for the Ising model in d = 2 and 
d = 3, respectively. Data cross around the known critical tem- 
peratures T c = 2.269 and T c = 4.511. In d = 2, 6 ~ 0.625 for 
T < T c and fl c ~ 0.74, while in d = 3, ~ 1.3 for T < T c and 
C ~ 1.4. In the inset of the upper panel the quantity X(t, T) 
is plotted over r for different temperatures. In the numerical 
simulations we consider systems ol N = L d spins, with L = 400 
for d = 2 and L = 128 for d = 3. 



diminishes by increasing r for T > T c , grows until it con- 
verges to Xq when T < T c , and is r independent at T c . In 
the upper panel of Fig. Q]the quantity X(t,T) is plotted 
over the temperature for different values of r, for the two- 
dimensional Ising model evolving via Glauber dynamics. 
One clearly observes that the curves intersect in a narrow 
region giving T c = 2.268 ± 0.002 in agreement with the 
analytical result. 

In the case of d = 3, since 8 C ~ 1.4 > 1 [IT)] , 
the integral J Q dt t~ 0a diverges, making X(t,T) use- 
less for extracting critical behaviors. One then can 
overcome this problem considering the "second moment" 

X 2 (t, T) = i Y^ftU^TJ > ^at * s ex P cc t G d to converge 
to (2 — # c )/(3 — C ) for the critical quench. In the lower 
panel of Fig. Q] the quantity ^(t, T) is plotted versus the 
temperature for different values of r. One again clearly ob- 
serves that the curves intersect in a narrow region giving 
T c = 4.510 ± 0.003, in agreement with previous numerical 
results. 



Fig. 2: (Color online). X(t,T) vs T for different values of r 
for the EA model in d = 2. In the inset T c (t;) is plotted vs n 
supporting the absence of a finite T c . 



As intermediate case between ferromagnets and spin 
glass systems, we have also considered the diluted Ising 
model. In this model the spin coupling Jij is chosen 
to be J with probability p and to be with probabil- 
ity 1 — p. Our results substantially agree with previous 
accurate estimates of T c in Ref. [7], for different choices 
of the parameter p. As a further test, we have consid- 
ered the d — 2 Edwards- Anderson model, where previous 
studies clearly indicate the absence of a phase transition 
at finite temperature {20Tl22j . More precisely, we have 
investigated systems with N = 400 2 spins and bimodal 
couplings Jjj = ±1 with equal probability. Results, plot- 
ted in fig 12 show that curves corresponding to different 
r do not intersect in the same point but the intersection 
points move towards the left by increasing r. More pre- 
cisely, in the inset of figHJ we plot the "critical tempera- 
ture" T c (Tj) identified from the intersection between the 
curves for X(n,T) and X(r,+i,T). We chose rj+i = an 
with Ti = 1000 and a = 1.35. The inset clearly shows that 
T c {Ti) decreases for increasing Ti, indicating an asymptotic 
convergence Tj — > at large times. This result supports 
the absence of a finite critical temperature. The same 
analysis performed for the other cases where a finite T c 
has been identified, gives a T c (rj) fluctuating around T c . 

Next we consider the case of the one-dimensional Ising 
spin glass with power- law decaying interactions |23j . The 
system is defined by the Hamiltonian H = — ■ JySjSj — 
h J^- Si, where the site i belongs to a ring of length L and 
h is a magnetic field. The sum is over all spins of the 
ring and = c(a)eij /rfj, where are chosen according 
to a Gaussian distribution with zero mean and standard 
deviation unity. The constant c(cr) is chosen to give a 
mean field transition temperature T^ IF = 1, namely 



1 = (T C MF ) 2 



E 



Jf, 



r 



(5) 
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Fig. 3: (Color online). The quantities Xi(t,T) and X(t,T) are plotted versus the temperature for different values of r, in 
the one-dimensional Ising spin glass with power-law interactions, for a — 0.55 and h — 0,0.1. In all cases the curves cross, 
indicating the presence of a finite critical temperature. 




Fig. 4: (Color online). The quantities X(t, T) is plotted over the temperature in the one-dimensional Ising spin glass with power- 
law interactions, for a — 0.55 and h = 0, 0.05, 0.1, 0.2. Curves correspond to different n obtained by the relation r»+i = 1.35t; 
and t\ — 1000. In the insets the plot T c (n) vs n gives a finite critical temperature for h = 0, 0.05, 0.1 and suggests T c = for 
h = 0.2. 
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Fig. 5: (Color online). Panel (a): the critical temperatures 
of the power law one-dimensional Ising spin glass measured in 
Refs. [26], [27] and in this study are plotted for several values 
of a in the absence of external field. Panel (b): the T — h phase 
diagram is reported for the same model with a = 0.75. 



where [. . .] av denotes an average over disorder. The dis- 
tance between two spins on the ring in terms of L is 
= (L/tt) sin(7r|z — j\/L). By varying the strength of 
the interaction through the parameter er, this model shows 
different behaviors [23j . In particular, for a = 0, taking 
c(ct) ~ 1/y/N, one recovers the Sherrington-Kirkpatrick 
model [24]. For a € [1/2,1] the system shows a fi- 
nite critical temperature, with a mean field-like region for 
a e [1/2, 2/3] and a non-mean field region for er e (2/3, 1]. 
In recent years, this model and its diluted version have 
been widely investigated in the literature [2"5l{2"9"] . focus- 
ing on the identification of a transition in presence of an 
external field, namely on the identification of the AT line. 
In these studies, T c for different values of a, both with 
and without the external field, has been measured with 
FSS analysis [25H29] , Contradictory results have been ob- 
tained in the non mean field-like region in presence of the 
external perturbation. When / 0, indeed, Leuzzi et 
al. |28j find a finite critical temperature, whereas no tran- 
sition has been observed by Katzgraber & Young [26l[29] . 
We turn to consider our results. In all our simulations we 
consider systems with N = 1024 spins and final times up 
to 10000 Monte Carlo steps. We have explicitly checked 
that no finite size effects are present. We first discuss re- 
sults for different values of a <G [0.5, 0.9] without external 
field. In particular, in the left panels of Figs. [3] and E] 
we plot X 2 (t, T) for a = 0.55, since 9 C = 1.09 ± 0.02 > 1, 
and X(t,T) for a = 0.75 versus T for different r. Curves 
clearly show intersection points, giving T c = 0.90 ± 0.02 
and T c = 0.61 ±0.02 for a = 0.55 and for a = 0.75, respec- 
tively. The critical temperature obtained from the same 
analysis for other values of a are reported in the left panel 
of Fig. [5] where previous results [26j[27] are also shown. 
We find that T c is a monotonously decreasing function of 
a, consistent with a linear decay Tc ~ 1 — 1.6(er — 0.5), 
for a G [0.5,0.8] and a faster decay for larger a. The 
monotonic decreasing behavior of T c with increasing a is 



expected, since larger values of a correspond to shorter in- 
teraction ranges. Conversely, the results of Ref. [55] show 
a decreasing linear behavior for a > 0.55, but are ex- 
pected to manifest a non-monotonic behavior approaching 
a = 0.5, where the mean field value T c — 1 is imposed by 
Eq. ([5]). We wish to notice that, for every choice of er, we 
always obtain a value of T c significantly smaller than the 
one obtained in Ref. [26]. The value of T c estimated in 
Ref. [57] are intermediate between our results and those 
of Ref. [13], but due to the large error bar, are compatible 
with both findings. 

Differences between the results of our method and those 
of Ref. [26] become more pronounced when ft, ^ 0. More 
precisely, in the mean field region, for a = 0.55 and h = 0.1 
we find T c = 0.71 ± 0.03, a value still smaller than that of 
Ref. [25] . T c = 0.96 ± 0.02. Conversely, an opposite trend 
is obtained in the non mean field region er > 2/3, where 
we always find a T c > 0, while no transition was obtained 
in Ref. [26]. Let us stress that we can consider a uniform 
field hi = h, whereas FSS analysis imposes the application 
of a spatially decorrelated (random) field. 

In particular, we focus on er = 0.75 and three values 
of h, h = 0.05,0.1,0.2. For h = and h = 0.1, we con- 
sider longer simulations up to 50000 Monte Carlo steps 
and N=2048 in order to avoid finite size effects. Results 
for X(T, t) are plotted in the main panels of Fig. [4] The 
behavior of T c (ri) vs n is plotted in the insets. Fig.[4]3 and 
Fig. 0J: show that curves intersect at a finite temperature 
indicating the existence of a phase transition in presence 
of an external field. We notice that the curves spreading at 
low temperature is less pronounced than for smaller a and 
other models. This can be attributed to the smaller value 
of 8 C . Indeed, from Eq. ^j, at fixed difference 9 C — 9 one 
has that Xo — X(t, T c ) is a decreasing function of 8 C . The 
measured values T c = 0.56±0.01 and T c = 0.43±0.01, for 
h = 0.05 and h — 0.1, respectively, are consistent with the 
expected trend of a decreasing T c for increasing h. The in- 
set of Fig.@Jl gives a non-constant T c (ri) that decreases at 
small Ti and tends to flattens only for the largest r^. The 
above trend suggests the absence of a finite T c for h = 0.02 
but does not exclude that T c (ri) asymptotically converges 
to T c G (0, 0.2) at large times. In the right panel of Fig. [S] 
we plot the T — h phase diagram showing the existence 
of the AT line, separating the spin glass phase from the 
paramagnetic one. 

In order to obtain more insights on the behavior of 
disordered systems, the same procedure can be carried 
out replacing C(t) in Eq.@ by a multi-spin correla- 
tion function. Nevertheless, the non linear susceptibil- 
ity Xi(t,t w ) = J2i,j( s i{ t ) s j{ t ) s i( t w)sj(t w )) usually consid- 
ered in the investigation of disordered systems |14j , is not 
properly suitable for this kind of study. Indeed, X4(M™) 
encodes the typical length scale L(t) of spatial correla- 
tion, and is expected to grow in time until a limit value 
that depends on £(T) or L{t w ) [3D]. More precisely, for 
L(t w ) < £(T),X(t,T) would be not affected by £(T) and, 
therefore, it would not be able to identify eventual diver- 
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gences of £(T). A way to overcome the above difficulty 
is to consider other multi-spin correlations, such as the 
second-order susceptibility considered in Ref. [31]. This 
quantity is intimately related to X4(t, t w ) and converges to 
an asymptotic value controlled by £(T), independently of 
t w . This quantity can be used in the described procedure. 
Difficulties, in this case, are related to huge fluctuations of 
this second order susceptibility that make the numerical 
evaluating very time-demanding. 

In conclusion, we have introduced an dimcnsionlcss 
quantity X(t,T) that is expected to become time inde- 
pendent at the critical temperature T c . This allows one 
to identify T c from the intersection of curves X(t,T) for 
fixed t and different T. The method has been tested in 
models where there exist accurate estimates of T c , as the 
Ising fcrromagnct and its diluted version. The study of 
the two-dimensional EA model confirms the absence of a 
transition at finite temperature. The method has been 
then applied to the one-dimensional Ising spin glass with 
power law decaying interactions for different choices of a 
and h. Results for ft, = give T c values always smaller than 
those obtained by static FFS methods in Ref. [2B], with 
differences that are larger for smaller a. In particular for 
a = 0.55 we obtain T c — 0.90 ± 0.02, a value smaller than 
the one of Ref. [26], where T c = 1.03 ±0.03, unexpectedly 
above the mean- field value T c = l. The study in presence 
of a finite perturbation h > indicates the existence of a 
finite critical temperature also in the non mean-ficld-likc 
region a > 2/3, in agreement with a replica symmetry 
breaking scenario. 
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